Scientific Machine Learning has become nowadays a hot topic of research, whose aim is to take profit from the strengths of the third (computational simulation) and fourth (data science) pillars of science, while alleviating their limitations. The blending of Machine Learning (ML) tools and physics-based Partial Differential Equations [1] or thermodynamics ideas [2] has become a competent approach for predicting the state of a system and gaining some insight about its inherent structure. However, the interpretability of the different model components is still obscure. Indeed, model knowledge is either restricted to some parameters whose value must be learned, to prescribed black box components of the whole model architecture, or to data structural invariants that do not necessarily coincide with physical variables. For instance, Physically-Guided Neural Networks with Internal Variables have been used for model selection [3] and to unravel the nonlinear nature of the material response [4], but there is still no decision rules based on statistical criteria as there are for conventional parametric fitting from structured data. To overcome this limitation, we introduce the so-called Reliability Aware Physically Guided Neural Network with Internal Variables (RAPGNNIV). This concept allows the assimilation of data uncertainty during the learning process, which is propagated to the explanatory components. Therefore, the overall framework acts as a knowledge distiller, as it allows not only to recover the structure of the eventual constitutive model, relating internal and measurable variables, but also to evaluate its uncertainty, thus enabling model identification and selection using statistical tools. Moreover, by prescribing a finite family of model candidates, it is possible to evaluate model plausibility or inadequacy, or to perform hyper-robust predictions [5]. We show that our framework can unravel the nonlinear nature of constitutive equations from unstructured datasets, generated with heterogeneous stimuli, thus leveraging the power of ML for equations discovery, while keeping track of the constitutive equation’s uncertainty. The method transforms the unstructured sensor measurements into pieces of information, encoded in the predictive block, able to recreate the system response for new conditions, and into elements of scientific knowledge, encoded in the explanatory block, enabling a deeper interpretation of its structure and properties.